In this lecture we will introduce linear models and least squares linear regression. We will explain the basic parametric model formulation and how to optimize it using linear algebra and geometry rather than gradient descent.
In this notebook we will be using the following packages:
import numpy as np
import pandas as pd
import plotly.offline as py
import plotly.express as px
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import plotly.figure_factory as ff
import cufflinks as cf
cf.set_config_file(offline=True, sharing=False, theme='ggplot');
import torch
from sklearn.linear_model import LinearRegression
# Support embedding YouTube Videos in Notebooks
from IPython.display import YouTubeVideo
This is the first of a series of lectures that are designed to be completed online rather than in class. This new lecture format is experimental and we ask for your feedback on Piazza.
The playlist for all the videos in this lecture can be found here:
Least Squares Linear Regression
We have designed this video notebook to be a series of short videos followed by simple questions (with answers provided) and demos. Note that all the solutions are provided and there are no grades for working through this notebook. However, we hope that by combining videos with activities this will actually be more engaging.
The following short video provides a recap of the modeling recipe (model, loss, optimize). If you are comfortable with this you can safely skip this video.
YouTubeVideo("u0YycwQ814c")
What are the three stages of modeling used in this class?
In the following video we introduce the multiple linear regression setting.
YouTubeVideo("21PN6a3s22I")
To motivate the multilinear regression setting we will work with the classic diamonds dataset.
data = pd.read_csv("data/diamonds.csv.zip", index_col=0)
data.head()
Each record in the dataset corresponds to a single diamond. The fields are
Fair, Good, Very Good, Premium, and Ideal)J (worst) to D (best).I1 (worstt), SI2, SI1, VS2, VS1, VVS2, VVS1, IF (best)We are interested in predicting the price of a diamond given it's characteristics. Therefore, we would like to fit a model:
$$ f_\theta\left(\text{Diamond's Characteristics}\right) \rightarrow \text{Price} $$Suppose for now we focus only on the diamond's carat, depth, and table values. In this case, we would want a model of the form:
$$ f_\theta\left(\textbf{carat}, \textbf{depth}, \textbf{table}\right) \rightarrow \textbf{price} $$We could express this as a linear model of the form:
$$ f_\theta\left(\textbf{carat}, \textbf{depth}, \textbf{table}\right) = \theta_0 + \theta_1 * \textbf{carat} + \theta_2 * \textbf{depth} + \theta_3 * \textbf{table} $$Where $\theta_0$ is the intercept parameter and $\theta_1$, ..., $\theta_3$ are the "slopes" with respect to each of the covariats.
Focusing on just the carat, depth, and table features and the price our dataset looks like:
data[['carat', 'depth', 'table', 'price']]
In lecture, we saw that the predictions for the entire data set $\hat{\mathbb{Y}}$ can be computed as:
$$ \hat{\mathbb{Y}} = \mathbb{X} \theta $$The covariate matrix $\mathbb{X} \in \mathbb{R}^{n,d+1}$ has $n$ rows corresponding to each record in the dataset and $d+1$ columsn corresponding to the origina $d$ columns in the data plus an additional 1s column vector. For example, the first five rows of $\mathbb{X}$ for this dataset would look like:
| 0 | 1 | 2 | 3 | |
|---|---|---|---|---|
| 0 | 1.0 | 0.23 | 61.5 | 55.0 |
| 1 | 1.0 | 0.21 | 59.8 | 61.0 |
| 2 | 1.0 | 0.23 | 56.9 | 65.0 |
| 3 | 1.0 | 0.29 | 62.4 | 58.0 |
| 4 | 1.0 | 0.31 | 63.3 | 58.0 |
We can extract our covariate matrix and add a column of all 1s:
n,d = data[['carat', 'depth', 'table']].shape
print("n = ", n)
print("d = ", d)
X = np.hstack([np.ones((n,1)), data[['carat', 'depth', 'table']].to_numpy()])
X
We can define the linear model in python using numpy. Here I will assume the input is in row vector form so this function will work with a single row or the entire matrix. Note in lecture we discuss the case where a single row is stored in column vector form (as this is the case in most math textbooks).
def linear_model(theta, xt):
return xt @ theta # The @ symbol is matrix multiply
If we guess a value for $\theta$ we can make a prediction using our model:
theta_guess = np.random.randn(d+1, 1)
theta_guess
Making a single prection for the 3rd row:
linear_model(theta_guess, X[2,:])
Making a prediction for all the rows at once:
linear_model(theta_guess, X)
Copying the notation directly from the video lecture:
Y_hat = X @ theta_guess
Y_hat
The next video defines and minimizes the least squares loss by using a geometric argument.
YouTubeVideo("zkJ3CULN8Wk")
The Least Squares Regression loss is simply the average squared loss.
$$ L(\theta) = \frac{1}{n}\sum_{i=1}^n \left( \mathbb{Y}_i - \mathbb{X}_i \theta \right)^2 $$Here $\mathbb{Y}_i$ is the $i$th observed response (diamond price) and $\mathbb{X}_i$ is the row vector coresponding to the $i$th row of the covariates with an added 1 at the beginning.
We need to define the response vector $\mathbb{Y}$
Y = data[['price']].to_numpy()
Y
Directly writing the average squared loss:
def squared_loss(theta):
return np.mean((Y - X @ theta)**2)
squared_loss(theta_guess)
Using matrix notation:
def squared_loss(theta):
return ((Y - X @ theta).T @ (Y - X @ theta)).item() / Y.shape[0]
squared_loss(theta_guess)
In the video lecture we defined the normal equation as:
$$ \mathbb{X}^T \mathbb{X} \hat{\theta} = \mathbb{X}^T \mathbb{Y} $$If we solve for $\hat{\theta}$ this is the minimizer of the squared loss with respect to our data.
If we assume $\mathbb{X}^T \mathbb{X}$ is invertible (full rank) then we get the following solution:
$$ \hat{\theta} = \left( \mathbb{X}^T \mathbb{X} \right)^{-1} \mathbb{X}^T \mathbb{Y} $$The following video provides a better intuition of the shape of each of the terms in the solution to the normal equation as well as how to better solve it.
YouTubeVideo("9vWR-19WQKU")
While it is not as numerically stable or efficient. We can compute $\hat{\theta}$ by direction using matrix inversion:
from numpy.linalg import inv, solve
theta_hat = inv(X.T @ X) @ X.T @ Y
theta_hat
theta_hat = solve(X.T @ X, X.T @ Y)
theta_hat
The solve can be implement slightly more efficiently by avoiding the inversion and subsequent matrix multiply. For matrices which are less well behaved, solve can also be more numerically stable.
We can compare the $\hat{\theta}$ to our earlier random guess. Notice here I will take the square root to make the loss in units of US dollars.
print("theta_hat", np.sqrt(squared_loss(theta_hat)))
print("theta_guess", np.sqrt(squared_loss(theta_guess)))
Better than random guess but not by much.
Now that we have estimated the model parameters $\hat{\theta}$ we can now also predict the price $\hat{\mathbb{Y}}$ for each of our diamonds
Y_hat = X @ theta_hat
In previous lectures we have plotted the residual. We can do that for each of the dimensions but in general it will be difficult to visualize the residual directly.
fig = px.scatter(x = X[:,1], y = Y - Y_hat,opacity=0.2)
fig.update_xaxes(title="Carat")
fig.update_yaxes(title="Residual")
fig
fig = px.scatter(x = X[:,2], y = Y - Y_hat,opacity=0.2)
fig.update_xaxes(title="Depth")
fig.update_yaxes(title="Residual")
fig
fig = px.scatter(x = X[:,3], y = Y - Y_hat,opacity=0.2)
fig.update_xaxes(title="Table")
fig.update_yaxes(title="Residual")
fig
A better way to visualize higher dimensional regression models is to compare the Predicted Value against the Observed value
fig = px.scatter(x = Y.flatten(), y = Y_hat.flatten(), opacity=0.2)
fig.add_trace(go.Scatter(x=[0,20000], y=[0, 20000], name="Y_hat=Y"))
fig.update_xaxes(title="Y")
fig.update_yaxes(title="Y_hat")
fig
What do we see? We are under estimating expensive diamonds! Perhaps we should consider the cut, color, and clarity.
In the next notebook lecture we will see how to improve this model by applying feature engineering techniques to incorporate categorical data.